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We present an exact result for the non-adiabatic transition probability and hence the defect density 
in the final state of a one-dimensional Kitaev model following a slow quench of the parameter J_, 
which estimates the anisotropy between the interactions, as J-(t) ~ — |t/r|. Here, time t goes 
from — oo to +oo and r defines the rate of change of the Hamiltonian. In other words, the spin 
chain initially prepared in its ground state is driven by changing J_ linearly in time up to the 
quantum critical point, which in the model considered here occurs at at t = 0, reversed and then 
gradually decreased to its initial value at the same rate. We have thoroughly compared the reverse 
quenching with its counterpart forward quenching i.e., J_ ~ t/r. Our exact calculation shows that 
the probability of excitations is zero for the wave vector at which the instantaneous energy gap is 
zero at the critical point J_ = as opposed to the maximum value of unity in the forward quenching. 
It is also shown that the defect density in the final state following a reverse quenching, we propose 
here, is nearly half of the defects generated in the forward quenching. We argue that the defects 
produced when the system reaches the quantum critical point gets redistributed in the wave vector 
space at the final time in case of reverse quenching whereas it keeps on increasing till the final time 
in the forward quenching. We study the entropy density and also the time evolution of the diagonal 
entropy density in the case of the reverse quenching and compare it with the forward case. 

PACS numbers: 73.43.Nq, 05.70.Jk, 64.60.Ht, 75.10.Jm 



I. INTRODUCTION 



Studies of quantum phase transitions in quantum 
many particle systems have always been a fascinating 
area of research in condensed matter physics P, ■ While 
a plethora of theoretical works have been performed on 
statics of quantum phase transitions (QPT), the dynam- 
ics of a quantum system passing through a quantum crit- 
ical point_has c_aught_ the attention of_ researchers only 
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2E 12J, |25|, |26|, |27[ . Understanding 
quantum dynamics happens to be a challenging problem 
as the physics of equilibrium quantum phase transitions 
gets coupled to the non-equilibrium dynamics of corre- 
lated systems. Theoretical techniques used in studying 
the above dynamics are often borrowed from similar stud- 
ies i n ^cj uantum optics, e.g., the Landau-Zcncr transition 

The dynamical evolution can be initiated in a quan- 
tum system either by a sudden change of a parameter 
in the Hamiltonian which is called a sudden quench 
or by a slow quenching of a parameter 0, ||. The ef- 
fect of the passage through a quantum critical point is 
manifested in the eventual dynamics of the system. The 
relaxation time of the system, defined as the time 
taken by the system to come back to its equilibrium state 
after a small perturbation, diverges at the critical point 
as £t ~ 8~ vz where 8 measures the deviation from the 
quantum critical point and v and z are the corrcspond- 
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ing correlation length and dynamical critical exponents, 
respectively. The divergence of the relaxation time is an 
artifact of vanishing energy gap A (~ £ t _1 ) between the 
ground and the first excited state of the Hamiltonian near 
a quantum critical point. This divergence of the relax- 
ation time forces the system to be infinitely sluggish near 
the critical point so that it takes effectively an infinite 
time to respond to any change in the external parame- 
ters thereby causing excitations. The recent discovery of 
ultra cold atoms which has facilitated the experimental 
implementation of various Hamiltonian models (30l . [3l| 
and thus the verification of results of quantum dynam- 
ics, has enormously accelerated the theoretical research 
in this field. Here, we are interested in a slow and linear 
variation of the quenching parameter and estimate var- 
ious quantities such as density of defects and the local 
entropy density in the final state of the system following 
a quench, as a function of the quenching rate r. 

It is well known that for a d-dimcnsional system which 
is initially prepared in its ground state and is quenched 
through a quantum critical point by linearly varying a 
parameter as t/r, the density of defects (n) satisfies 
the Kibble-Zurek(KZ) scaling || H, El, H ^ given by 
n ~ T -dv/(vz+l) wnere v anc i z are the critical exponents 
defined above. The Kibble-Zurek scaling has been veri- 
fied in various exactly solvable spin models 0, H H, EH EH 
and in a system of interacting bosons undergoing super- 
fluid to insulator transitions @. The above KZ scal- 
ing relations gets modified when the system is quenched 
through a multicritical point |25| . across a gapless phase 
[ltl EH , along a gapless line [2ll.l23l] or for quenching with 
a non-linear rate [17J. Studies on quenching dynamics 
have also been generalized to quantum spin chains with 
quenched disorder [Io| , to systems in presence of white 



2 



noise 14 1, to systems with infinite range interactions [22| . 
to an open system coupled to a heat bath [l9||, and also 
to quantum spin chains driven by an oscillatory magnetic 
field [26| . The effect of edge states on the defect produc- 
tion has also been studied [27l ]. It is worth mentioning 
here that the defect production has been studied experi- 
mentally for a rapid quench in a spin-1 bose condensate 

In this paper, we study the effect of the reversal of 
the quenching path right at the quantum critical point 
on the density of defects. This is achieved by increas- 
ing the quenching parameter from time — oo to its value 
at the quantum critical point and the bringing it back 
at the same rate to its initial value at the final time 
i.e., t — > oo. We call this quenching scheme as reverse 
quenching whereas the other scheme in which the quench- 
ing parameter is monotonically increased from time — oo 
to +00 through a quantum critical point will be referred 
to as the forward quenching scheme. We would also oc- 
casionally use the term half quenching for the case when 
the quenching is stopped at the quantum critical point. 

At the outset, let us discuss the motivation behind our 
study. The forward quenching scheme has been applied 
extensively for the entire range in time from t = — 00 
to 00 and defects generated in the final state has been 
estimated. In the present work, we drive the system lin- 
early right up to the quantum critical point and then let 
it retrace its path. In our calculation, we find that the 
defect generated up to time t = is approximately half 
of the defect generated for a full forward quenching. - 
We also seek answer to the question that how the defects 
generated in the first half of the quenching get altered 
under reversal of the path! We address questions like 
do we have a similar scaling form for the defect or how 
does the magnitude of the defects in the final state, i.e. 
at t — ► 00 change under reverse quenching?. In the pro- 
cess, we also provide an exact solution of the Schrodinfer 
equation to find the non-adiabatic transition probability 
for the reverse quenching. 

We also compare our work with that in ref. [2(| where a 
quantum XY spin chain is repeatedly swept through the 
quantum critical points by varying the magnetic field at 
a linear rate between —00 and 00 such that the reversal 
of path takes place far away from the critical points. In 
the present work, the parameter is increased only upto 
the quantum critical point where it is reversed to trace 
back its path. To employ the reverse quenching scheme 
in an appropriate way, we study the dynamics of a one- 
dimensional Kitaev model where the quantum critical 
point occurs at t = 0. Secondly, in ref. [20j, the defect 
density after each repetition is estimated using a recur- 
sive relation for the non-adiabatic transition probabili- 
ties where the (rapidly varying) cross terms are ignored 
because they vanish upon integration over the wave vec- 
tors. On the other hand, we present here an exact so- 
lution of the transition probabilities which include in- 
terference terms although the results arc in fairly good 
agreement at least qualitatively with the one cycle case 



discussed in ref [20( It is also to be noted that the reverse 
quenching dynamics has been studied for a generic two 
level system in ref. [35| and the excitation probability 



has been calculated within the framework of perturba- 
tion theory. On the other hand, we here generalize the 
quenching scheme to a many particle system and solve 
the Schrodingcr equations exactly. 

The paper is organized in the following way: the model 
and the quenching scheme are discussed in section II. The 
main results for the non-adiabatic transition probability, 
defect density and the entropy density are presented in 
section III. We summarize our results in the concluding 
section whereas the calculational details are provided in 
the appendix. 



II. THE MODEL AND THE QUENCHING 
SCHEME 



Two-dimensional Kitaev model defined on a honey- 
comb lattice described by the Hamiltonian [36[ 

H= ]T + J20 V n _ u o v nl + J 3 < i < i+1 )(l) 

n+!-even 

where n and I defines the column and row indices of 
the lattice has been studied extensively due to its exact 
solvability by Jordan- Wigncr transformation [37] . The 
rich phase diagram of this model has a gapless phase, 
through which the quenching dynamics has been studied 
recently [Tr|. The one dimensional version of the Kitaev 
model (with J3 = 0) given by the Hamiltonian [36|, H|| 



H 



N 

£ 

n=l 



( ~Y X '3l> t y u \ 



(2) 



where n refers to the site index, exhibits a quantum phase 
transition at J\ = J 2 . The above Hamiltonian @, which 
is the model of interest in this work, can be exactly diag- 
onalized by standard Jordan Wigner transformation [371 ] 
as defined below 



2n+l- 



(3) 



Here a n and b n are independent Majorana fcrmions at 
site n [lfj. They satisfy the relations like 



{b m , b n } = 25 m>n , {a m , b n } = 0. 



(4) 



Substituting for af. L and in terms of Majorana fermions 
followed by a fourier transformation, Hamiltonian ([2]) can 
be written as 

7T 

H = 2iJ2[b{a k (Ji + J2e lk )+alb k (-J 1 - J 2 e~ tk )} (5) 

k=0 
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where the fourier component a,k, satisfying the stan- 
dard anticommutation relations {ak,a\.i} — &k.k' an d 
{afe, afc'} = 0, is defined as 



afce 



ikn 



-ikn] 



k=0 



+\l jj[ao + oj + a*(-l) n + 4(-l) r 



(0) 



The sum over k in Eq. ((6]) goes only for half the Brillouin 
zone as a' n s are Majorana fermions. By defining ipk = 
(a,k,bk), the Hamiltonian ([5]) can then be rewritten in a 
simpler form as 



below: In the limit of large t (t — > ±oo). the eigenvectors 
of the Hamiltonian are 



eifc - 



e 2 fc 



lk k 

v /2(l + sin( fc /2)) [COS 2 lfe) + (1 + Sin 2)'^] 



where \e±k) is the ground state in the limit t — > — oo. A 
unitary transformation generated by the matrix U con- 
structed from the above eigenvectors leads to the final the 
Hamiltonian suitable for the present form of quenching 
and is given by 



(7) 



where 



-Ji - J 2 e" ifc 
J i + J 2 e 4fc 



(8) 



The above Hamiltonian can be diagonalized where the 
eigenvalues arc given by 



= ±2 
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2 Ji J 2 cos k. 



Clearly, the gap vanishes at J\ = ±J 2 for k = it and 
respectively with the critical exponents v and z both 
being equal to unity. Feng, Zhang and Xiang [3^ showed 
that this vanishing energy gap signals a topological phase 
transition between the two phases of the model at J± < 
J 2 and Ji > J 2 • Interestingly, this model can be mapped 
to a one-dimensional transverse Ising model by a duality 
transformation [H, S3, • 

In terms of a new set of basis functions given by 



ipik = -j= ( \ I and ip 2 k 
the above Hamiltonian can be recast to the form 



1 



1 

—i 



H, 



J++J- 



J+-J- 



cos k 



J+-J- 



sin k 



J+-J- 



sin k 



J++J- 



J+-.J- 



cos k 



H'k 



U^H k U 

J-(t) sin(fc') 
J + cos(fc') 



J + cos(fc') 
-J_(t) sin(fc') 



(10) 



where the time dependence is now entirely shifted to di- 
agonal terms of the Hamiltonian. The quantum critical 
point is at J_ = for the mode k' = 7r/2. Also, the mode 
k' in Eq. (p~0|) is half of mode k in Eq {8}. Henceforth, we 
will refer fc' as /c and appropriately redefine the Brillouin 
zone. The presence of a single QCP precisely at t = 
renders the analytical calculation easier and so we chose 
Kitaev model over other exactly solved models for the 
present study. We note that the results of this model can 
be extended to any other Jordan- Wigncr solvable models. 



where J± = J\ ± J 2 . We study the dynamics of the 
spin chain by varying the term J_ of Hamiltonian ([9|) 
using the quenching rule J_ = — 1-| where t varies from 
—00 to +00. Here, the quantum critical point occurs 
at t = and it is at this point where the parameter 
J_ is reversed to bring it back to its initial value at the 
final time. It is to be noted that the off-diagonal terms 
in the Hamiltonian ([9]) becomes time-dependent in the 
present quenching scheme making the analytical solution 
difficult. However, the situation can be easily saved by 
making an appropriate unitary transformation as shown 



III. RESULTS 

In this section, we shall present the main results of 
this work. The 2x2 reduced Hamiltonian matrix given 
in Eq. (1101) can be interpreted as a Landau Zcner Hamil- 
tonian |28l I29I ] where the diagonal elements arc the two 
bare (diabatic) energy levels which approach each other 
and the off-diagonal element is the minimum gap be- 
tween the instantaneous levels of the Hamiltonian. At 
time t = 0, the energy gap between the instantaneous 
energy levels vanish for the mode k = tt/2 signaling a 
(g\ quantum phase transition mentioned above. We shall as- 
sume that the system is prepared in its initial ground 
state \e\k) at t — ► —00. At any instant t during the evo- 
lution, a general state vector \ipf~(t)) can be expressed as 
\ipk{t)) = ci,k(t)|eifc) + c 2 ,fc(i)|e 2 fc) where Ci tk (t) (i=l,2) 
denote the time-dependent probability amplitude for the 
bare state \eik)- 

The Schrodingcr equation describing the evolution of 
the system is 

d 

igl^.kit) = 2J_sin(fc)ci ifc (t) + 2cos(fc)c 2:fe (i) 
d 

«^c 2 , fc (i) = -2J_sin(fc)c 2 ,fc(i)+2cos(fc)ci, fc (i), (11) 
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with initial conditions c\.k{— oo) = 1, C2,fc(— oo) = and 
we have set J + = 1. The nonadiabatic transition prob- 
ability in the final state is given by |c2,fc(+oo)| 2 . The 
above Schrodinger equations are solved exactly and the 
probability of excitation for the fc-th mode, pk, is 



Pk (t -> oo) 



i' 1 



r(l-ia/2) .r(l/2-ia/2) 



r(l + ia/2) r(l/2 + ia/2) 



(12) 



where a = rcos 2 (fc)/sin(fc). The density of defects can 
be obtained by integrating the probability of excitations 
Pk over the Brillouin zone and is given by 



i r i r 

7T / Pk(t ->■ oo) dk = - / 



p k {t^X,)dk. (13) 



The parameter a measures the effective rate of driving. 
It is a, not r which determines the diabatic (a — > 0) 
and adiabatic (a — > oo) limit [1, H2|- The gap varies 
with wave vector k, so does a and for the modes close 
to the critical mode (k = it/2), a ~ k 2 r. We defer the 
calculational details to the Appendix. 

Let us first analyze the exact expression given in 
Eq. (|12p in different limits and compare it with pk ob- 
tained by direct numerical integration of the Schrodinger 
equation (fTTj) . It should be noted that the first expression 
in the modulus squared term of Eq. (|12D is 



r(l-za/2) _ T(z) _ T(z) 
T(l + ia/2) ~f(I) ~ f^J 

which is a unit vector with argument — 29i where 0\ is the 
argument of T(l + ia/2). Similarly the second expression 
in the modulus squared term is also a unit vector with 
angle — 2# 2 and hence modulus squared term in Eq. (CO 
reduces to 



2 + 2sin(26» 2 -26)i) 



(14) 



Therefore, the final expression for the probability of 
excitations is 



p k (t -> oo) 



1 



(1 



-2-na 



) x |2 + 2sin(20 2 - 20i)|(15) 



Using the properties of V function [431] . we have 



^2^ + (^2(TT7)- tan % 



2(1 + 0' 



• t_XJ 

y ?l 

f-^ 2(1/2 



, -tan _1 (— r) 116) 

(1/2 + V 2(l/2 + /) ; / ' 



where 4>( z ) 1S the well known Digamma function [43j|. The 
analytical expression for p k is obtained by substituting 
8 1 and 02 in Eq (fT5)) . The result obtained by numerical 
integration of the Schrodinger equation is presented in 
Fig. ([T]) where we also plot the analytical expression after 
substituting numerically obtained values of 0\ and 02 in 

Eq. m 




FIG. 1: Variation of pk vs k for r = 1. The data with '+' 
sign represent the numerical solution whereas the dashed line 
corresponds to the analytical expression given in Eq. 1121 As 
explained in the text, the region near k = n/2 varies linearly 
with a — rcos 2 fc/sinfc whereas that away from k = ty/2 as 
1/a 2 . The inset also shows the intersection of the two limits 
at a particular mode ko- The dotted line in the inset goes as 
7ra and the thin line falls as l/16a 2 . 



The behavior of pk as a function of the wave vector k 
can be explained by making resort to the Landau-Zcncr 
interpretation discussed before. For the modes close to 
k = (which are away from the critical mode k = 7r/2), 
the minimum gap is relatively large, or more precisely 
Af.r >> 1. Hence, these modes evolve adiabatically re- 
maining close to the instantaneous ground state through- 
out the quenching. For the critical mode k = it/2, gap 
is zero or in other words the relaxation time (inverse of 
the gap) is diverging which results to the complete freez- 
ing of dynamics. The system stays in its initial ground 
state throughout the evolution which also happens to be 
the ground state at t — > +oo for the present scheme of 
quenching. We therefore encounter a situation where the 
mode for which instantaneous energy gap is zero has si- 
multaneously zero excitation probability which is in con- 
trast to the forward case where the probability of exci- 
tations is unity for the critical mode. Similarly, for the 
modes near k = tt/2 where the gap is still very small, the 
system stays closer to the initial state due to large relax- 
ation time leading to a final state similar to the ground 
state. We therefore conclude that the transition proba- 
bility vanishes in either limits k — > and k — > tt/2 and a 
peak is expected at a wave vector k Q lying somewhere in 
the middle as shown in Fig. (1). 

The instantaneous excitation at an instant t is defined 
as the probability of finding the system in the instanta- 
neous excited eigenstatc of the Hamiltonian (10) . The 
variation of the instantaneous excitation probability as a 
function of time presented in Fig.® reflects the expla- 
nation presented above. 

In the limit of small a (i.e., either r small or k — > tt/2), 
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t 

FIG. 2: Instantaneous excitation probabilities vs time for 
three different modes. The mode k = 1.49 is closer to the 
7r/2 mode and shows a diabatic behavior i.e., the small gap 
and large relaxation time makes the system unable to change 
appreciably from the initial state and therefore the instanta- 
neous excitation decreases in magnitude for t > when it is 
retracing its path. On the other hand, the system tries to fol- 
low the instantaneous ground state for k = 1.3 which is in the 
adiabatic limit and the excitation keeps on increasing till the 
effect of finite gap persists. Finally, for the mode k = 1.55 
which is closest to the critical mode k — n/2, decrease in 
instantaneous excitation for t > is prominently shown. 

wc get a simplified expression 

26> 2 - 26»i = -2a In 2. (17) 

Substituting Eq. (fTT]) in Eq (fT5]l . we get the expression 
for pk in the small a limit as 

p k = -(1 - e~ 27ra ) x (2 - 2sin(2aln2)). (18) 

which correctly predicts the curve for small a along with 
the peak of the curve. 

It would be useful to calculate pk in the two extreme 
limits, namely, a — > and a — > oo. 

The a — ► limit can be obtained directly from Eq. ( 
I12p whereas the large a limit is obtained from using the 
asymptotic expansion of the Gamma function in Eq (|12[) . 
Thus, we have 

Pk ~ 7rcv for a — > 

Pk ~ -—^ for a — > oo (19) 

which is also shown in the inset of Fig. 1. As discussed 
already, the a — > behavior is expected near the critical 
mode k — ir/2 whereas large a behavior can be observed 
for the modes k away from 7r/2. It is also interesting 
to note that the exact solution we present here reduces 
to the adiabatic transition probability scaling 1/a 2 in 
the limit of large a as expected from the quantum me- 
chanical adiabatic theorem [2^, [44| ■ This feature is more 
transparent in Fig. [3] where the variation of pk with r 
for two different values of k, one near k = tt/2 and the 
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FIG. 3: Variation of logp^ vs logr for two different k values. 
Fig (a) corresponds to k=1.3 where 1/t 2 behavior is expected. 
The dots are the numerically obtained values where as the 
fitted line has a slope -2. Fig (b) (inset), on the other hand 
is for k = 1.56 where increase linearly with r. Once again 
a log-log plot shows a slope of 1 as expected from theory. 

other away from k = 7r/2, is shown depicting the linear 
increase with r and decrease as 1 /r 2 respectively. 

The mode ko at which the peak in pk occurs is approx- 
imately given by the point of intersection of two limiting 
behaviors given in Eq. 15 and is given by 

cos 2 (fco) = , 1 q /3 l 
sin(fco) 167T t 

The value of a at ko is given by 

a(k = k ) = i/l/16ir 

which is independent of r. This implies that the maxi- 
mum value of pk is independent of r as shown in FiglH A 
rough estimate of this value can be obtained from Eq (fT8|) 
after substituting a at ko and is found to be close to 0.25. 

Let us shift our attention to estimating the density of 
defects n in the final state at t — > oo. The variation 
of density of defects with the rate of quenching r is ob- 
tained by integrating the probability of excitations over 
the Brillouin zone given by Eq. (13). We find that n(r) 
as a function of r shows a peak at a particular quench- 
ing rate To and eventually follows a 1/ y/r decay for very 
large r. The 1/y/r behavior is justified by noting the 
fact that for very large r, ko shifts towards 7r/2 where 
cos 2 (fc)/ sin(fc) ~ k 2 . In this large r-limit therefore, pk 
scales as pk ~ Pk(Tk 2 ) resulting to the Kibblc-Zurck scal- 
ing of the defect density given as n ~ 1/y/r. Fig. 5 shows 
the variation of density of defects with r in the large r 
limit with a 1/y/r behavior whereas the inset of Fig. 5 
corresponds to the n vs r behavior for the entire range 
of t depicting the peak as described above. 

An interesting observation is that for the present 
quenching scheme, the density of defects in the final state 
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FIG. 4: The variation of pk for different values of r showing 
that the maximum value of pk is independent of r 
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FIG. 5: The main part of the figure shows the variation of 
density of defects n for relatively higher values r for reverse 
quenching. Also plotted is half times the density of defects 
produced while forward quenching and it is clear that reverse 
quenching is close to half of the linear quenching. Inset shows 
n vs. r for a wider range of r where a peak is observed for a 
relatively smaller value of r 
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FIG. 6: a) A comparison between probability of excitations 
for reverse quenching, forward quenching and pk(t — 0). is 
presented The range of k is appropriately chosen for the clar- 
ity of presentation. Fig. 6b)shows that the defect generated 
in the reverse case, which is the area under the pk vs. k curve, 
is approximately equal to the defect generated at t = 0. To 
highlight this, we have doubled pt for the reverse case and 
overlapped one of its peaks with pk{t = 0) by appropriate 
shifting of the x-axis. 



is close to half of that in the forward quenching i.e., the 
case where </_ is linearly quenched from — oo to oo, see 
Fig This is because the maximum value of pk in the 
reverse case is one-fourth that of the forward case making 
the area under one of the peaks to be close to one fourth 
that of the linear. 

It is also illustrative to compare the non-adiabatic 
transition probability Pk(t — > +oo) as a function of k 
for reverse quenching, forward quenching as well as half 
quenching i.e., pk(t = 0). Fig. [6ja) suggests that the den- 
sity of defects (area under the pk vs k curve) does not 
change appreciably in reverse quenching as compared to 
the half quenching. However, there is a reorganization 
of pk in the wave- vector space keeping the density of de- 
fect nearly constant. To justify the above statement, we 



have doubled the peaks in the reverse case and appropri- 
ately shifted the x-axis to match one of its peaks with 
Pk (t — > 0) in Fig EJb) . The peak is found to match al- 
most identically to Pk(0). In the passing, we note from 
Fig. (6a), that for very large r when only the modes 
close to the critical mode contributes to the defects, the 
density of defects in the forward case is double that of 
half quenching. 

We now shift our focus to the local von-Neuman en- 
tropy density @ defined by 

S=-- f (p k lnp fe + (1 - p k ) Hi - p k )) dk (20) 

generated in the reverse quenching process. The entropy 
density is small for small as well as large values of r 
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FIG. 7: Comparison of entropy of reverse and forward quench- 
ing. 



and attains the maximum value at a characteristic time- 
scale. A comparison between the entropy generated in 
the forward (Sp) and the reverse case (Sr) is shown in 
Fig. [7J Although qualitatively the curves are similar, 
an interesting difference is to be pointed out: Sr is less 
than Sf in the limit of small r whereas Sr exceeds Sf 
for large r. This observation can be explained as follows: 
the integrand in Eq (|20p is maximum when pk = 0.5. In 
reverse case, pk always remains less than 0.5, i.e., it never 
reaches the maximally disordered state. For the forward 
case in the small t (non-adiabatic) limit, pk is close to 
0.5 in a larger region and hence entropy is large. In other 
words, the final state is more locally ordered following a 
reverse quenching than a forward quenching in the limit 
r — > 0. On the other hand, for the large r limit in the 
forward case, pk increases sharply near the critical mode 
k = tt/2 and non-negligible only for wave vectors close to 
7r/2 resulting to relatively smaller entropy density. 

In a recent work, Barankov and Polkovnikov [45[ have 
proposed the concept of diagonal entropy given by 

Sd = — ^2 P nn ^ n p™ 1 

n 

where p nn is the n-th diagonal clement of the density 
matrix describing the system. One can interpolate it to 
obtain a time dependent diagonal entropy where p nn (t) 
is the diagonal elements of the density matrix in the in- 
stantaneous eigen basis. In our case, pn(t) = 1 — Pk(t) 
and P22{t) = Pk(t) where the excitations for each mode k 
are calculated in the instantaneous eigen basis. Here, we 
compare the evolution of the diagonal entropy in the re- 
verse and forward cases. We find that in the forward case 
the diagonal entropy increases monotonically with time 
and eventually saturates to the asymptotic value corre- 
sponding to the von-Neuman entropy whereas in the re- 
verse case, a dip is observed immediately after the critical 
point. 

The experimental realization of the Kitaev model has 
been proposed recently in systems of ultracold atoms and 
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time 

FIG. 8: Variation of diagonal entropy with time for r = 10. 



molecules trapped in optical lattices [3 II ] - In this pro- 
posal, each of the couplings can be independently tuned 
using different microwave radiations. Once this is estab- 
lished, one can experimentally verify the reverse quench- 
ing case by looking at the defect density which actually 
corresponds to the number of bosons in the wrong spin 
state. It is also possible to investigate the spatial corre- 
lation function of the operator ib n a n + r where a n and b n 
are Majorana fermions as defined before [l6| . This spa- 
tial correlation function depends on pk which we have 
already obtained for the reverse case. Then the evo- 
lution of defect correlations can be detected by spatial 
noise correlation measurements as discussed in Ref. 46. 
Also, qualitative testing of reverse quenching can be done 
by varying a magnetic field in spin gap dimer compounds 
such as BaCuSi20g which undergo a singlet-triplet quan- 
tum phase transition at a critical field B c . Our results 
suggest that the density of defects in the reverse case 
which correspond to residual singlets obtained by mag- 
netization measurement would be close to half of the for- 
ward case. 



IV. CONCLUSIONS 

In conclusion, we study the dynamics of a one- 
dimensional Kitaev model using a reverse quenching 
scheme, by increasing the anisotropy parameter J_ lin- 
early up to its quantum critical point at t = following 
which J_ is decreased at the same rate to bring it back 
to its initial value. We provide an exact solution of the 
Schrodingcr equation and estimate the density of defects 
and the local entropy density in the final state. Com- 
parison of the reverse quenching results to those of the 
corresponding forward quenching case leads to a few in- 
teresting observations. In the reverse case, the Landau- 
Zener transition probability pk vanishes for the critical 
mode k = tt/2 at which the instantaneous energy gap 
vanishes at t = whereas Pk is maximum for the same 
mode in the forward case [16| . We show that pk increases 
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linearly with a in small a (diabatic) limit whereas our re- 
sult retrieves the expected 1/a 2 fall of p^ in the large a 
(adiabatic) limit as predicted from quantum mechanical 
adiabatic theorem; a is the effective rate parameter as 
defined in the text. Interestingly, the density of defects 
in the reverse case is close to half that of forward quench- 
ing. We have also compared the half quenching case with 
reverse quenching. A close inspection of Pk as a function 
of k as shown in Fig. 6 suggests that in the reverse case, 
Pk gets reorganized in wave vector space with a peak at 
k = fco 7^ 7r/2 keeping the density of defects approxi- 
mately same as half quenching case. The value of this kg 
shifts to ir/2 for t — > oo. We also show that the max- 
imum value of the transition probability is independent 
of t. The local entropy density of the final state and 
time evolution of the diagonal entropy density are also 
explored. The possibility of experimental realization of 
the reverse quenching scheme using atoms trapped on 
optical lattices has also been pointed out. 
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APPENDIX A: DETAILS OF EXACT 
CALCULATION 

We present here an outline of the calculational details 
leading to an exact solution for the reverse quenching 
case generalizing earlier studies d, [H, H3| to the present 
case. Let us first consider a general Hamiltonian in the 
basis |1) and |2) as shown below 



H 



ex A 

A* e 2 



where ex and e 2 are the two bare energy levels (diagonal 
elements) varying as ~ t/2r and —t/2r respectively (see 
Eq. {TO) for comparison). If \ij}(t)) = c x (i)|l) + c 2 (t)\2), 
then one can write down the Schrodinger equation for 
C\(t) and c 2 (t). But before that we redefine c\{t) and 
c 2 (t) as follows 



ci(t) = £ 1 (t)e-' i/ -- e ' 1<Jt ' 

ca(t) = c 2 {t)e- iS -^ dt ' . 

The Schrodinger equation for c 2 is 

i^-ut) = Acim e -^-»(^*')-«(*')) dt '. 

ot 

One more transformation of the form 



(Al) 



(A2) 



(A3) 



helps us to write the equation for c 2 in terms otU 2 in the 
following form 



(A4) 



where we have substituted e\ — e 2 = t/r. Now redefining 
a new variable 



t 



-iir/A 



one gets, 

d 2 1 



-)U 2 {z) = 



(A5) 



where m = iA 2 T. By all these transformations, we are 
able to recast the Schrodinger equation in the form of We- 
ber Differential equation [47| whose solutions are linear 
combination of well known Weber functions .D_ m _i(iz) 
and £)_ TO _i(— iz) i.e., 

U 2 (z) = aD_ m _ x {iz) + bD-m-xi-iz) (A6) 
or going back to the notation of cx(t) and c 2 (t), 

\m) = ^[d t -^][aD_ m _x(iz) + bD_ m _x(-iz))\l) 
+ [aD_ m _x{iz) + bD_ m _ x {-iz)}\2). (A7) 

But the initial condition demands that at t — > — oo, 
!?/>(£)) ~ |1) forcing U 2 {z) to be a function of only 
D— m —x(— iz) as D- m -x(—iz) goes to zero at t — » — oo 
but D- m -x(iz) does not as can be seen from the follow- 
ing asymptotic form of Weber functions 



D n {z)~e-^z n 



2 71 " n-Rir \z 2 -n-1 

T(-n) 6 6 

irr 5irr 
for — < arg(z) < — — 
4 w 4 



where r is either 1 or —1, and 
D n (z) — e~i* 2 z n for \arg(z)\ < ^ 



(A8) 



(A9) 



Therefore, a — and the form of b can be obtained by 
these asymptotic forms along with the initial condition, 
which gives b = Ay / re^7 A r . 

Hence, Eq (|A7|) . after using the derivative of Weber 
function, is 

IV (*<0)) = ri A!T e«T[(m+l)fl_ m _ 2 H Z ) 
- iz£>_ m _!(-iz)]|l) + A^e-^ A2r D_ m _x{-iz)\2) 

At t = 0, the wavefunction is 



m = 0)) = e-t A2 V^ ^ 2 



re 4 



At 



r(l/2 + m/2) 

7T 2- m ' 2 



I) 



2 r(l + m/2) 



|2) (A10) 
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which is obtained by using the following property 

lim D m {s) = 2 m '* n .f" j- + 0(a). 
s-+o r(l/2 - m/2) 

The wave function at t > 0, which has the effect of re- 
versing, should match with the wave function for t < 
at t = 0. In the reverse case, the parameters m and z for 
i > are redefined as 



-iA 2 T and z' 



-it 
7? 



(All) 



With these redefined k' and z', starting from Eq (|A7[) . 
the wavefunction for i < are matched with that of t > 
at t = to obtain coefficients a and b 



1 



Av^e" 



rA-T<- 



r(l-m/2) r(l/2-m/2) 



r(l + m/2) n/2 + m/2 



6= -A % /7e"7 A T 2~ m x 



r(l-m/2) .r(l/2-m/2)" 

~T~ % ~ 



r(l + m/2) IT/2 + m/2 



We know that at t — > oo, the excited state is |2) and hence 
the coefficient \c2\ 2 defines the excitation probability and 
is equal to 



c 2 (t —> oo) ~ lim aL>_ m /_i(iz') + 6D_ m /_i(-iV) (AI2) 



Once again, using the expression for b and asymptotic 
expansion of Weber function in the definition of |c2 (t — > 
oo)| 2 with properly identifying r and A for a Kitaev 
model, we get Eq (|T^|) . 
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